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Overview of Multirate Signal Processing 


Digital transformation of the sampling rate of signals, or signal processing 
with different sampling rates in the system. 


Applications 


1. Sampling-rate conversion CD to DAT format change, for example. 

2. Improved D/A, A/D conversion oversampling converters; which 
reduce performance requirements on anti-aliasing or reconstruction 
filters 

3. FDM channel modulation and processing bandwidth of individual 
channels is much less than the overall bandwidth 

4. Subband coding of speech and images Eyes and ears are not as 
sensitive to errors in higher frequency bands, so many coding schemes 
split signals into different frequency bands and quantize higher- 
frequency bands with much less precision. 


Outline of Multirate DSP material 


1. General Rate-changing System 

2. Integer-factor Interpolation and Decimation and Rational-factor Rate 
Changing 

. Efficient Multirate Filter Structures 

. Optimal Filter Design for Multirate Systems 

. Multi-stage Multirate Systems 

. Oversampling D/As 

. Perfect-Reconstruction Filter Banks and Quadrature Mirror Filters 
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General Rate-Changing Procedure 


This procedure is motivated by an analog-based method: one conceptually 
simple method to change the sampling rate is to simply convert a digital 
signal to an analog signal and resample it! ([link]) 
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The problems with this scheme are: 


1. A/D, D/A, filters cost money 
2. imperfections in these devices introduce errors 


Digital implementation of rate-changing according to this formula has three 
problems: 


1. Infinite sum: The solution is to truncate. Consider sinc (t) ~ 0 for 
t < ti, t > to: Then mT, — nTo < ty and mT) — n7To > te which 
implies 


N2 
fn) = xo(n) sincr: (mT, — nTp) 
n=N, 


Note: This is essentially lowpass filter design using a boxcar window: 
other finite-length filter design methods could be used for this. 


. Lack of causality: The solution is to delay by max {||} samples. 
The mathematics of the analog portions of this system can be 
implemented digitally. 

Equation: 
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Equation: 
sin ( (mT, —nT)) 
x1(m) = De ro(n =n) T'Smax{Ty,1 4} 
= SO% .x0(n) sincr (mT; — nTo) 


So we have an all-digital formula for exact digital-to-digital rate 
changing! 

. Cost of computing sincz (mT — nT): The solution is to 
precompute the table of sinc (¢) values. However, if 7 is nota 
rational fraction, an infinite number of samples will be needed, so 
some approximation will have to be tolerated. 


Note: Rate transformation of any rate to any other rate can be 
accomplished digitally with arbitrary precision (if some delay is 
acceptable). This method is used in practice in many cases. We will 
examine a number of special cases and computational improvements, 
but in some sense everything that follows are details; the above idea is 
the central idea in multirate signal processing. 


Useful references for the traditional material (everything except PRFBs) are 
Crochiere and Rabiner, 1981 and Crochiere and Rabiner, 1983. A more 
recent tutorial is Vaidyanathan; see also Rioul and Vetterli. References to 
most of the original papers can be found in these tutorials. 


Interpolation, Decimation, and Rate Changing by Integer Fractions 


Interpolation: by an integer factor L 


Interpolation means increasing the sampling rate, or filling in in-between samples. Equivalent to 
sampling a bandlimited analog signal L times faster. For the ideal interpolator, 
Equation: 


Xo(Lw) if |u| < = 
Oif F< |wl<a 


Xi(w)={ 


We wish to accomplish this digitally. Consider [link] and [link]. 
Equation: 
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denoted as 
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The DTFT of y(m) is 
Equation: 


Y(w) = Dm=vy(mje 
= Vint to(mje WP") 


eae 2(m)en ven) 


= Xo (wL) 


Since Xq(w’) is periodic with a period of 2x, Xo(Lw) = Y(w) is periodic with a period of 2* (see 
Uink]). 
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By inserting zero samples between the samples of x9(n), we obtain a signal with a scaled frequency 
response that simply replicates Xo(w’) L times over a 27 interval! 


Obviously, the desired x;(m) can be obtained simply by lowpass filtering y(m) to remove the 


replicas. 
Equation: 
x1(m) = y(m)*hr(m) 
Given 
lif jwol)< 
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In practice, a finite-length lowpass filter is designed using any of the methods studied so far ({link]). 
Interpolator Block Diagram 


rae eee 


Decimation: sampling rate reduction (by an integer factor M) 


Let y(m) = ¢o(Lm) ([link]) 
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That is, keep only every Lth sample ([link]) 
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In frequency (DTFT): 
Equation: 
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Now DTFT [5° 6(n — Mk)] = 20 So X(k)6(w, — 22%) for |w| < 7 as shown in homework 
#1, where X(k) is the DFT of one period of the periodic sequence. In this case, X(k) = 1 for 
M-— us 
k € {0,1,...,M —1} and DTFT [> d(n — Mk)] = 20 opty 6(w, — 7). 
Equation: 
DTFT [zo(n)]*DTFT [> 5(n — Mk)] = Xo(w")*2x HG" 5 (w, — 288) 
= 2 Jt Kol!) (2m Dk" 6 (wr = wr — 2BE)) a 
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so Y(w) = 34g X0(- — 22#) ie., we get digital aliasing .((link]) 
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Usually, we prefer not to have aliasing, so the downsampler is preceded by a lowpass filter to 
remove all frequency components above |w| < +7 ([link)). 


Rate-Changing by a Rational Fraction L/M 


This is easily accomplished by interpolating by a factor of Z, then decimating by a factor of MZ 
({link]). 
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The two lowpass filters can be combined into one LP filter with the lower cutoff, 


1 if |u| <4, 
3 -| max{L,M} 


0 if mabey Sw <7 


Obviously, the computational complexity and simplicity of implementation will depend on 1: 2/3 
will be easier to implement than 1061/1060! 


Efficient Multirate Filter Structures 
Rate-changing appears expensive computationally, since for both 


decimation and interpolation the lowpass filter is implemented at the higher 
rate. However, this is not necessary. 


Interpolation 


For the interpolator, most of the samples in the upsampled signal are zero, 
and thus require no computation. ({link]) 
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Form = L | = | +m mod L and p =m mod L, 
Equation: 


zi(m) = So™? vy. Arp(m)y(m) 
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Pictorially, this can be represented as in [link]. 


y(m) or 
Bh 
T rate L'T rate 
0 0 


These are called polyphase structures, and the g,(7) are called polyphase 
filters. 


Computational cost 
If h(m) is a length-N filter: 
N IN computations 


¢ No simplification: Th — 


_1_ \ computations 
iy sec 


e Polyphase structure: (z+ _ 7 where L is the 


number of filters, - is the taps/filter, and a is the rate. 


Thus we save a factor of L by not being dumb. 


Note: For a given precision, N is proportional to L, (why?), so the 
computational cost does increase with the interpolation rate. 


Note: Can similar computational savings be obtained with IIR structures? 


Efficient Decimation Structures 


We only want every Mth output, so we compute only the outputs of 
interest. ([link]) 
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Polyphase Decimation Structure 
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The decimation structures are flow-graph reversals of the interpolation 
structure. Although direct implementation of the full filter for every Mth 
sample is obvious and straightforward, these polyphase structures give 
some idea as to how one might evenly partition the computation over M 
cycles. 


Efficient L/M rate changers 


Interpolate by Z and decimate by MM ((link)). 
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Combine the lowpass filters ({link]). 
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We can couple the lowpass filter either to the interpolator or the decimator 
to implement it efficiently ({link]). 
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Of course we only compute the polyphase filter output selected by the 
decimator. 
Computational Cost 


Every T; = “tT seconds, compute one polyphase filter of length 2, or 


+ + N multiplies 


Ti “1 ~ MT, second 


However, note that NV is proportional to max {L, M}. 


Filter Design for Multirate Systems 


The filter design techniques learned earlier can be applied to the design of filters in multirate 
systems, with a few twists. 


Example: 

Design a factor-of-Z interpolator for use in a CD player, we might wish that the out-of-band error 
be below the least significant bit, or 96dB down, and < 0.05 % error in the passband, so these 
specifications could be used for optimal LD, filter design. 


In a CD player, the sampling rate is 44.1kHz, corresponding to a Nyquist frequency of 22.05kHz, 
but the sampled signal is bandlimited to 20kHz. This leaves a small transition band, from 20kHz to 
24.1kHz. However, note that in any case where the signal spectrum is zero over some band, this 
introduces other zero bands in the scaled, replicated spectrum ([link]). 
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So we need only control the filter response in the stopbands over the frequency regions with nonzero 
energy. ([link]) 
Y(w) 
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The extra "don't care" bands allow a given set of specifications to be satisfied with a shorter-length 
filter. 


Direct polyphase filter design 


Note that in an integer-factor interpolator, each set of output samples z;(Ln + p), 

p = {0,1,..., L — 1}, is created by a different polyphase filter g,(n), which has no interaction 
with the other polyphase filters except in that they each interpolate the same signal. We can thus 
treat the design of each polyphase filter independently, as an 4 Jength filter design problem. 
({link]) 


Each g,(m) produces samples x;(Ln + p) = x (n + 2), where n + + is not an integer. That is, 
Gp() is to produce an output signal (at a Tp rate) that is x(n) time-advanced by a non-integer 


advance +. 


The desired response of this polyphase filter is thus 
Hp,p(w) = er 
for |w| < 7, an all-pass filter with a linear, non-integer, phase. Each polyphase filter can be designed 


independently to approximate this response according to any of the design criteria developed so far. 
Exercise: 


Problem: What should the polyphase filter for p = 0 be? 


Solution: 


A delta function: ho(n) = 6(n’) 


Example: 
Least-squares Polyphase Filter Design 


¢ Deterministic x(n) Minimize 


Given x(n) = x(n)*h(n) and za(n) = z(n)*ha(n). Using Parseval's theorem, this becomes 
Equation: 


min {72 (l2(n) — ea(n)|)?} 


min { 2 J”, (|X(w)H(w) — X(w)Ha(w)|)? dw} 


min { i %, |H(w) — Ha(w)|(|X(w)))? dw} 


This is simply weighted least squares design, with (|X(w)|)” as the weighting function. 
e stochastic X(@) 
Equation: 


min { E|(|e(n) — 2a(n)|)°| } 


B | (|e(n)*((n) — ha(n)) 1)? 
= min {2 f%, (|Ha(w) — H(w)|)?Sz0(w) dw} 


Sz (w) is the power spectral density of x. 


Sza(w) = DTFT [ree(*)| 
Pe (k) = Bla(k a )a(0)| 
Again, a weighted least squares filter design problem. 
Exercise: 
Problem: Is it feasible to use IIR polyphase filters? 


Solution: 


The recursive feedback of previous outputs means that portions of each IIR polyphase filter 
must be computed for every input sample; this usually makes IIR filters more expensive than 
FIR implementations. 


Multistage Multirate Systems 


Multistage multirate systems are often more efficient. Suppose one wishes 
to decimate a signal by an integer factor M, where M is a composite 
integer M = M,Mp2...M, = ty M,. A decimator can be implemented 
in a multistage fashion by first decimating by a factor Mj , then decimating 
this signal by a factor Mo , etc. ({link]) 

Multistage decimator 
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Multistage implementations are of significant practical interest only if they 
offer significant computational savings. In fact, they often do! 


The computational cost of a single-stage interpolator is: 


N_ taps 
MT sec 


The computational cost of a multistage interpolator is: 


ee 
MiTo MMT MT 


The first term is the most significant, since the rate is highest. Since 

N; « M,; for a lowpass filter, it is not immediately clear that a multistage 
system should require less computation. However, the multistage structure 
relaxes the requirements on the filters, which reduces their length and 
makes the overall computation less. 


Filter design for Multi-stage Structures 


Ostensibly, the first-stage filter must be a lowpass filter with a cutoff at ae 
to prevent aliasing after the downsampler. However, note that aliasing 
outside the final overall passband |w| < +7 is of no concern, since it will 
be removed by later stages. We only need prevent aliasing into the band 

|w| < +7; thus we need only specify the passband over the interval 


|w| < 3, and the stopband over the intervals w € aa Ss geet ae | 


M, M? M, 
fork € {1,...,M — 1}. ([link)) 
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Of course, we don't want gain in the transition bands, since this would need 
to be suppressed later, but otherwise we don't care about the response in 
those regions. Since the transition bands are so large, the required filter 
turns out to be quite short. The final stage has no "don't care" regions; 
however, it is operating at a low rate, so it is relatively unimportant if the 
final filter turns out to be rather long! 


L-infinity Tolerances on the Pass and Stopbands 


The overall response is a cascade of multiple filters, so the worst-case 
overall passband deviation, assuming all the peaks just happen to line up, is 


1 +4, =]]1+&, 
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So one could choose all filters to have equal specifications and require for 
each-stage filter. For 6,,, <1, 
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Alternatively, one can design later stages (at lower computation rates) to 
compensate for the passband ripple in earlier stages to achieve 
exceptionally accurate passband response. 


6, remains essentially unchanged, since the worst-case scenario is for the 
error to alias into the passband and undergo no further suppression in 
subsequent stages. 


Interpolation 


Interpolation is the flow-graph reversal of the multi-stage decimator. The 
first stage has a cutoff at + ((link]): 
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However, all subsequent stages have large bands without signal energy, due 
to the earlier stages ([link]): 
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The order of the filters is reversed, but otherwise the filters are identical to 
the decimator filters. 


Efficient Narrowband Lowpass Filtering 


A very narrow lowpass filter requires a very long FIR filter to achieve 
reasonable resolution in the frequency response. However, were the input 
sampled at a lower rate, the cutoff frequency would be correspondingly 
higher, and the filter could be shorter! 


The transition band is also broader, which helps as well. Thus, [link] can be 
implemented as [link]. 
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and in practice the inner lowpass filter can be coupled to the decimator or 
interpolator filters. If the decimator and interpolator are implemented as 
multistage structures, the overall algorithm can be dramatically more 
efficient than direct implementation! 


DFT-Based Filterbanks 


One common application of multirate processing arises in multirate, multi- 
channel filter banks ([link]). 


x(n) 


One application is separating frequency-division-multiplexed channels. If 
the filters are narrowband, the output channels can be decimated without 
significant aliasing. 


Such structures are especially attractive when they can be implemented 
efficiently. For example, if the filters are simply frequency modulated (by 

~(:%#4n) 
e€ 


FFTs! 


) versions of each other, they can be efficiently implemented using 


Furthermore, there are classes of filters called perfect reconstruction 
filters which are of finite length but from which the signal can be 
reconstructed exactly (using all M channels), even though the output of 
each channel experiences aliasing in the decimation step. These types of 
filterbanks have received a lot of research attention, culminating in wavelet 
theory and techniques. 


Uniform DFT Filter Banks 


Suppose we wish to split a digital input signal into NV frequency bands, 
uniformly spaced at center frequencies w, = 2th forO<k< N—-1. 
1 if jul < + 
. Bandpass 
0 otherwise 


filters can be constructed which have the frequency response 


Consider also a lowpass filter h(n), H(w) ~ { 


from 


The output of the kth bandpass filter is simply (assume h(n) are FIR) 
Equation: 


; 2nkm ) 


a(n)*hy(n) = soo} a(n —m)h(m)e— 
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This looks suspiciously like a DFT, except that IZ 4 N, in general. 
However, if we fix M = N, then we can compute all y;,(7) outputs 
simultaneously using an FFT of x(n — m)h(m): The 

kth FFT frequency output = yx(n)! So the cost of computing all of 
these filter banks outputs is O[N log NJ, rather than NV 2\ per a given n. 
This is very useful for efficient implementation of transmultiplexors 
(FDM to TDM). 

Exercise: 


Problem: 
How would we implement this efficiently if we wanted to decimate the 


individual channels y;,(n) by a factor of N, to their approximate 
Nyquist bandwidth? 


Solution: 


Simply step by NV time samples between FFTs. 
Exercise: 
Problem: 


Do you expect significant aliasing? If so, how do you propose to 
combat it? Efficiently? 


Solution: 
Aliasing should be expected. There are two ways to reduce it: 


1. Decimate by less ("oversample" the individual channels) such as 
decimating by a factor of x. This is efficiently done by time- 
stepping by the appropriate factor. 

2. Design better (and thus longer) filters, say of length DN. These 
can be efficiently computed by producing only N (every Lth) 
FFT outputs using simplified FFTs. 


Exercise: 


Problem: 


How might one convert from N input channels into an FDM signal 
efficiently? ([Link]) 


Note: Such systems are used throughout the telephone system, 
satellite communication links, etc. 


Solution: 


Use an FFT and an inverse FFT for the modulation (TDM to FDM) 
and demodulation (FDM to TDM), respectively. 


Quadrature Mirror Filterbanks (QMF) 


Although the DFT filterbanks are widely used, there is a problem with aliasing in the decimated channels. 
At first glance, one might think that this is an insurmountable problem and must simply be accepted. 
Clearly, with FIR filters and maximal decimation, aliasing will occur. However, a simple example will show 
that it is possible to exactly cancel out aliasing under certain conditions!!! 


Consider the following trivial filterbank system, with two channels. ({link]) 


Note @(n) = x(n) with no error whatsoever, although clearly aliasing occurs in both channels! Note that 
the overall data rate is still the Nyquist rate, so there are clearly enough degrees of freedom available to 
reconstruct the data, if the filterbank is designed carefully. However, this isn't splitting the data into separate 
frequency bands, so one questions whether something other than this trivial example could work. 


Let's consider a general two-channel filterbank, and try to determine conditions under which aliasing can be 
cancelled, and the signal can be reconstructed perfectly ({link]). 


Let's derive Z(n), using z-transforms, in terms of the components of this system. Recall ([link]) is 
equivalent to 


x(n)-»[H@)}» y(n) 


and note that ({link]) is equivalent to 
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and ([link]) is equivalent to 
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Y (z) is derived in the downsampler as follows: 


Yizj= » z(Mm)z"™ 
Let n = Mm andm = 7,7, then 
¥(z) = a x(n) » d(n — Mp)z ™ 
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Armed with these results, let's determine X(z) <= 2(n). ((link]) 
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and 


La(2) = 5 Fi(2)Hi(2)X(2) + 5 Fi(2)Ea(—2)X(—2) = 5 Fle) Hal) X(2) + 5 Fi (2)Ba(—2) X(-2) 


Finally then, 
Equation: 


X(2z) = Ua(z) + La(z) 
= ¥ (Ho(z)Fo(z)X(z) + Ho(—2)Fo(2)X(—z) + Hi (2)Fi(z)X(z) + Ai(-z) Fi(z)X(—2)) 
= 3 (Ho(2)Fo(z) + Hi(2)Fi(2))X(z) + 3 (Ho(—2)Fo(z) + Hi(—2)Fi(z))X(—z) 


Note that the X(—z) —> X(w + 7) corresponds to the aliasing terms! 
There are four things we would like to have: 


1. No aliasing distortion 

2. No phase distortion (overall linear phase — simple time delay) 
3. No amplitude distortion 

4. FIR filters 


No aliasing distortion 


By insisting that Hy)(—z) F(z) + H(—z)F,(z) = 0, the X(—z) component of X(z) can be removed, and 
all aliasing will be eliminated! There may be many choices for Ho, H1, Fo, F) that eliminate aliasing, but 
most research has focused on the choice 


Fo(z) = Hi(—z) : Fi(z) = —Ho(—z) 
We will consider only this choice in the following discussion. 


Phase distortion 

The transfer function of the filter bank, with aliasing cancelled, becomes 

T(z) = + (Ao(z)Fo(z) + Ai(z)Fi(2)), which with the above choice becomes 

T(z) = + (Ho(z)Hi(—z) — Hi(z)Ho(—z)). We would like T(z) to correspond to a linear-phase filter to 
eliminate phase distortion: Call 


P(z) = H(z) Hi (—z) 
Note that 


1 
T(z) = 5 (P(2) — P(-2)) 
Note that P(—z) <= (—1)"p(n), and that if p(n) is a linear-phase filter, (—1)"p(n) is also (perhaps of the 
opposite symmetry). Also note that the sum of two linear-phase filters of the same symmetry (i.e., length of 
p(n) must be odd) is also linear phase, so if p(n) is an odd-length linear-phase filter, there will be no phase 
distortion. Also note that 


2p(n) if nis odd 
O if nis even 
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means p(n) = 0, when n is even. If we choose ho(n) and hi (7) to be linear phase, p(7) will also be linear 
phase. Thus by choosing ho(7) and hi (7) to be FIR linear phase, we eliminate phase distortion and get FIR 
filters as well (condition 4). 


Amplitude distortion 
Assuming aliasing cancellation and elimination of phase distortion, we might also desire no amplitude 
distortion ( |T'(w)| = 1). All of these conditions require 


T(z) = ; AH) -MoniDse 


where c is some constant and D is a linear phase delay. c = 1 for |T’(w)| = 1. It can be shown by 
considering that the following can be satisfied! 


2p(z) = 2cd(n — D) if nis odd 


E. = P = P = = 2 ae 
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Thus we require 
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Any factorization of a P(z) of this form, P(z) = A(z)B(z) can lead to a Perfect Reconstruction filter bank 
of the form 


H(z) = A(z) 
HH 1 (—z) = B(z) 
[This result is attributed to Vetterli.] A well-known special case (Smith and Barnwell) 


Ay(z) =—- (2-0?}4H9(-z")) 


Design techniques exist for optimally choosing the coefficients of these filters, under all of these constraints. 
Equation: 
Quadrature Mirror Filters 


Ay (z) = Ho(—z) & Ai(w) = Ho(m + w) = Hy (a —w) 
for real-valued filters. The frequency response is "mirrored" around w = +. This choice leads to 
T(z) = Ho*(z) — Ho?(—z): it can be shown that this can be a perfect reconstruction system only if 
Ho(z) = coz) + cyz7@™) 


which isn't a very flexible choice of filters, and not a very good lowpass! The Smith and Barnwell approach 
is more commonly used today. 


M-Channel Filter Banks 


The theory of M-band QMFBs and PRFBs has been investigated recently. 
Some results are available. 


Tree-structured filter banks 


Once we have a two-band PRFB, we can continue to split the subbands 
with similar systems! ([link]) 
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Thus we can recursively decompose a signal into 2. bands, each sampled at 
2 th the rate of the original signal, and reconstruct exactly! Due to the tree 
structure, this can be quite efficient, and in fact close to the efficiency of an 
FFT filter bank, which does not have perfect reconstruction. 


Wavelet decomposition 


We need not split both the upper-frequency and lower-frequency bands 
identically. ([link]) 
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low low-high frequencies 
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This is good for image coding, because the energy tends to be distributed 
such that after a wavelet decomposition, each band has roughly equal 
energy. 


